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Abstract: In this manuscript we provide an outline of the numerical meth¬ 
ods used in implementing the density constrained time-dependent Hartree-Fock 
(DC-TDHF) method and provide a few examples of its application to nuclear 
fusion. In this approach, dynamic microscopic calculations are carried out on a 
three-dimensional lattice and there are no adjustable parameters, the only input 
is the Skyrme effective NN interaction. After a review of the DC-TDHF theory 
and the numerical methods, we present results for heavy-ion potentials V{R), 
coordinate-dependent mass parameters M{R), and precompound excitation en¬ 
ergies E*{R) for a variety of heavy-ion reactions. Using fusion barrier penetra¬ 
bilities, we calculate total fusion cross sections cr(Fc,m.) for reactions between 
both stable and neutron-rich nuclei. We also determine capture cross sections 
for hot fusion reactions leading to the formation of superheavy elements. 
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1 Introduction 

The theoretical investigation of heavy-ion interaction potentials V (i?) is of fundamental 
importance for the study of fusion reactions at energies in the vicinity of the Coulomb 
barrier. Radioactive ion beam facilities enable us to study fusion reactions with exotic 
neutron-rich nuclei [T]. An important goal of these experiments is to study the effects of 
neutron excess {N — Z) on fusion. Another theoretical challenge is to provide guidance to 
fusion experiments leading to the formation of new superheavy elements. And finally, fusion 
of very neutron rich nuclei plays a major role in nuclear astrophysics where it determines 
the composition and heating of the crust of accreting neutron stars. A recent review of the 
state of heavy-ion fusion including the DC-TDHF method can be found in [5]. 

In the absence of a practical quantum many-body theory of barrier tunneling, sub-barrier 
fusion calculations are commonly performed in two steps; the calculation of a heavy-ion 
interaction potential V(R) and a subsequent calculation of tunneling through the barrier. 
Theoretical studies of fusion cross sections are often done by employing phenomenological 
methods such as the coupled-channels approach 0 ID in which one uses empirical ion- 
ion potentials or double-folding potentials calculated with “frozen density” approximation 
and low-lying excitations are incorporated using macroscopic models. These calculations 
ignore dynamical effects such as neck formation, particle exchange, pre-equilibrium GDR 
and other modes during the nuclear overlap phase of the reaction. During this phase of 
the collision the primary underlying mechanism is the dynamical change in the density 
along the fusion path which modifies the potential energy. Obviously, this density change 
is not instantaneous. For instance, it was shown in Ref. [S] that the development of a 
neck due to couplings to octupole phonons in '^^Ca+'^^Ga could take approximately one 
zeptosecond. As a consequence, the dynamical change of the density is most significant 
at low energy (near the barrier-top) where the colliding partners spend enough time in the 
vicinity of each other with little relative kinetic energy. At high energies, however, the nuclei 

* To be published as a commemorative ebook focusing on the field of nuclear reaction theory and on the 
work related to Prof. Joachim Maruhn. 



overcome the barrier essentially in their ground-state density. This energy dependence of 
the effect of the couplings on the density evolution was clearly shown in TDHF calculations 
by Washiyama and Lacroix [^. This naturally translates into an energy dependence of the 
nucleus-nucleus potential Clll], similar to what was introduced phenomenologically in the 
Sao-Paulo potential [9] . Consequently, the barrier corresponding to near barrier-top energies 
includes dynamical couplings effects and can be referred to as a dynamic-adiabatic barrier, 
while at high energy the nucleus-nucleus interaction is determined by a sudden potential 
which can be calculated assuming frozen ground-state densities. 

We have developed a dynamic microscopic approach, the density constrained time- 
dependent Hartree-Fock (DC-TDHF) method [10], to calculate heavy-ion interaction po¬ 
tentials V{R), mass parameters M{R), and precompound excitation energies E*{R) [TT] 
directly from microscopic TDHF dynamics. The idea of constraining the dynamical den¬ 
sity during a TDHF collision to find the collective path corresponding to the dynamical 
changes of the nuclear density was first introduced in a collaborative work of Prof. Maruhn 
in 1985 [T2|. This novel approach has facilitated the present calculations of sub-barrier fu¬ 
sion cross sections, capture cross sections for superheavy element production, and nuclear 
astrophysics applications. The theory has been applied to calculate fusion cross-sections 
for the fusion of light to heavy systems ElllMI]. In this paper we will summarize the 
DC-TDHF method and outline the numerical algorithms involved in solving the static and 
dynamic equations of motion. For a variety of heavy-ion reactions, we discuss numerical 
calculations of fusion / capture cross-sections and compare these to experimental data if 
available. 


2 Theoretical formalism 

The theoretical formalism for the microscopic description of complex many-body quan¬ 
tum systems and the understanding of the nuclear interactions that result in self-bound, 
composite nuclei possessing the observed properties are the underlying challenges for study¬ 
ing low energy nuclear physics. The Hartree-Fock approximation and its time-dependent 
generalization, the time-dependent Hartree-Fock theory, has provided a possible means to 
study the diverse phenomena observed in low energy nuclear physics including collective 
vibrations [22H25] and nuclear reactions [261 [27] . 

2.1 Time-dependent Hartree-Fock (TDHF) method 

Given a many-body Hamiltonian containing two and three-body interactions 

N N N 

H = ^ ^ ti -|- ^ ^ Vij -\- ^ ^ Vijk , ( 1 ) 

i i<j i<j<k 

the time-dependent action S can be constructed as 

S = f " dt < - ihdtmt) > . (2) 

Jti 

Here, $ denotes the time-dependent, many-body wavefunction, <i)(ri, r 2 ,..., ta; t), and ti 
is the one-body kinetic energy operator. General variation of S recovers the time-dependent 
Schrodinger equation. In TDHF approximation the many-body wavefunction is replaced by 
a single Slater determinant and this form is preserved at all times. The determinental form 
guarantees the antisymmetry required by the Pauli principle for a system of fermions. In 
this limit, the variation of the action yields the most probable time-dependent path between 
points ti and t 2 in the multi-dimensional space-time phase space 

<5S' = 0^$(t) = l>o(t) . (3) 

In practice is chosen to be a Slater determinant 

^o(i) = -^det\(l)x{r ,t)\ , 


(4) 



where ,t) are the single-particle states with quantum numbers A. If the variation in 
Eq.(^ is performed with respect to the single-particle states (j)*\ we obtain a set of coupled, 
nonlinear, self-consistent initial value equations for the single-particle states 

h{{4>f^})(t)x = X = l,...,N . (5) 

These are the fully microscopic time-dependent Hartree-Fock equations which preserve the 
major conservation laws such as the particle number, total energy, total angular momentum, 
etc. As we see from Eq.([^, each single-particle state evolves in the mean-field generated by 
the concerted action of all the other single-particle states. Static equations can be obtained 
from Eq.([^ by taking out a trivial phase from the single-particle states 

K{x^l})xx = eAXA 

• ( 6 ) 

In TDHF, the initial nuclei are calculated using the static Hartree-Fock (HE) theory. 
The resulting Slater determinants for each nucleus comprise the larger Slater determinant 
describing the colliding system during the TDHF evolution. Nuclei are assumed to move on 
a pure Coulomb trajectory until the initial separation between the nuclear centers used in 
TDHF evolution. Using the Coulomb trajectory we compute the relative kinetic energy at 
this separation and the associated translational momenta for each nucleus. The nuclei are 
then boosted by multiplying the HF states with 

—)■ exp(zkj • R)$j , (7) 

where is the HF state for nucleus j and R is the corresponding center of mass coordinate 

R=ix;r.. (8) 

^ i—1 

The Galilean invariance of the TDHF equations assures the evolution of the system without 
spreading and the conservation of the total energy for the system. In TDHF, the many-body 
state remains a Slater determinant at all times. The final state is a filled determinant, even 
in the case of two well separated fragments. This phenomenon is commonly known as the 
“cross-channel coupling” and indicates that it is not possible to identify the well separated 
fragments as distinct nuclei since each single particle state will have components distributed 
everywhere in the numerical box |28j . In this sense it is only possible to extract inclusive 
(averaged over all states) information from these calculations. 

2.2 DC-TDHF method 

The concept of using density as a constraint for calculating collective states from TDHF 
time-evolution was first introduced in Ref. [12) . and used in calculating collective energy 
surfaces in connection with nuclear molecular resonances in Ref. [29) . 

In this approach we assume that a collective state is characterized only by density p, and 
current j. This state can be constructed by solving the static Hartree-Fock equations 

^ 0 J (9) 

subject to constraints on density and current 

< ®pjl/3(r)l^pj > = P{r,t) 

< ^pjliWI^pJ > = j(r,t). 

Choosing p{Y,t) and j(r, t) to be the instantaneous TDHF density and current results in the 
lowest energy collective state corresponding to the instantaneous TDHF state |$(t) >, with 
the corresponding energy 


Ecoii{p{t),iit)) —< ®pj|R|'^’p,j > • 


(10) 



This collective energy differs from the conserved TDHF energy only by the amount of internal 
excitation present in the TDHF state, namely 


E*{t) = Etdhf — Ecoii{t) ■ ( 11 ) 

However, in practical calculations the constraint on the current is difficult to implement but 
we can define instead a static adiabatic collective state |$p > subject to the constraints 

< $p|/5(r)|$p > = p(r,t) 

< ^p|j(i-)I^P > = 0. 

In terms of this state one can write the collective energy as 

Ecoii = Ek^■n{p{t),}{t)) + Edc{p{Y: t)) , ( 12 ) 

where the density constrained energy Edc^ and the collective kinetic energy Ekin are defined 
as 


Edc 

Ekin 


< ^p\H\^p > 


where the index q is the isospin index for neutrons and protons {q = ri,p). 


2.3 DC-TDHF applications: ion-ion potential, mass parameter, 
and precompound excitation energy 

Recently, we have developed a new method to extract ion-ion interaction potentials di¬ 
rectly from the TDHF time-evolution of the nuclear system [10]. In the DC-TDHF approach 
the TDHF time-evolution takes place with no restrictions. At certain times t or, equivalently, 
at certain internuclear distances R{t) the instantaneous TDHF density is used to perform a 
static DFT energy minimization while constraining the proton and neutron densities to be 
equal to the instantaneous TDHF densities. This means we allow the single-particle wave 
functions to rearrange themselves in such a way that the total energy is minimized, subject 
to the TDHF density constraint. In a typical DC-TDHF run, we utilize a few thousand time 
steps, and the density constraint is applied every 10 — 20 time steps. 

In essence, DC-TDHF provides us with the TDHF dynamical path in relation to the 
multi-dimensional static energy surface of the combined nuclear system. In this approach 
there is no need to introduce constraining operators which assume that the collective motion 
is confined to the constrained phase space. In short, we have a self-organizing system which 
selects its evolutionary path by itself following the microscopic dynamics. 

We refer to the minimized energy as the “density constrained energy” Ey)c ■ From Eq. 
is is clear that the density constrained energy plays the role of a collective potential, except 
for the fact that it contains the binding energies of the two colliding nuclei. One can thus 
define the ion-ion potential as (TOj 

V{R{t)) = Eucipir, t)) - Ea, - Ea, , (13) 

where Ea^ and Ea 2 are the binding energies of two nuclei obtained from a static Hartree- 
Fock calculation with the same effective interaction. For describing a collision of two nuclei 
one can label the above potential with ion-ion separation distance R(t) obtained during 
the TDHF time-evolution. This ion-ion potential V{R) is asymptotically correct since at 
large initial separations it exactly reproduces Vcouiomb(,Rmax)- All of the dynamical features 
included in TDHF are naturally included in the DC-TDHF potentials. These effects include 
neck formation, particle exchange , internal excitations, and deformation effects to all 

order, among others. Couplings between relative motion and internal structures are known 
to affect fusion barriers by dynamically modifying the densities of the colliding nuclei. The 


effect is expected to be stronger at energies near the barrier top, where changes in density 
have longer time to develop than at higher energies. This gives rise to an energy dependence 
of the barriers as predicted by modern time-dependent Hartree-Fock (TDHF) calculations [S]- 

m- 

The TDHF evolution also provides us with a coordinate-dependent mass M{R) which 
may be obtained using energy conservation for a central collision 


M{R) = 


2[F;c.m. - vm 
R^ 


(14) 


where the collective velocity R is directly obtained from the TDHF evolution. The in¬ 
dependence of this mass at lower energies is very similar to the one found in constrained 
Hartree-Fock calculations [33] with a constraint on the quadrupole moment. On the other 
hand, at higher energies the coordinate dependent mass essentially becomes flat, which is 
again a sign that most dynamical effects are contained at lower energies. The peak at small 
R values is due to the fact that the center-of-mass energy is above the barrier and the 


denominator of Eq. (14) becomes small due to the slowdown of the ions. 

The fusion barrier penetrabilities Tl^Ec.^.) can be obtained by numerical integration 
of the Schrodinger equation for the collective distance coordinate R, using the heavy-ion 
potential V{R) with coordinate dependent mass parameter M{R). Alternatively, we can 
instead use the constant reduced mass fi and transfer the coordinate-dependence of the 
mass to a scaled potential V (R) using the well known exact point transformation [13 [33] 




(15) 


The potential V(R), which includes the coordinate-dependent mass effects differs from the 
V{R) only in the interior region of the barrier. Further details can be found in Ref. |17j . 

The fusion barrier penetrabilities Tl(i?c.ni.) are obtained by numerical integration of the 
Schrodinger equation 


_ 2^J, dR^ 


L{L + l)fi^ 


+ V{R) - Ec.ni. 


HR) = 0 , 


(16) 


using the incoming wave boundary condition (IWBC) method |34j . IWBC assumes that 
once the minimum of the potential is reached fusion will occur. In practice, the Schrodinger 
equation is integrated from the potential minimum, Rmin, where only an incoming wave 
is assumed, to a large asymptotic distance, where it is matched to incoming and outgoing 
Coulomb wavefunctions. The barrier penetration factor, T]^{Ec.m.)j is the ratio of the incom¬ 
ing flux at i?niin to the incoming Coulomb flux at large distance. Here, we implement the 
IWBC method exactly as it is formulated for the coupled-channel code CCFULL described 
in Ref. |3S]. This gives us a consistent way for calculating cross-sections at energies below 
and above the barrier via 


a/(Ac.m.) = ^ ■ (17) 

L=0 


At energies above the barrier either the DC-TDHF method or direct TDHF calculations 
can be used to determine the fusion cross-sections. The precompound excitation energy 


E*{R{t)) can be computed from Eq. II 


The DC-TDHF calculations mentioned above utilize the potential V{R), which includes 
the coordinate-dependent mass effects, to find the cross-sections using Eq. Recently, a 
systematic study was performed to calculate the potential and coordinate dependent mass 
for non-central TDHF collisions [5]. 


2.4 Skyrme interaction 

Almost all TDHF calculations have been done using the Skyrme energy density func¬ 
tional. The Skyrme energy density functional contains terms which depend on the nuclear 










density, p, kinetic-energy density, r, spin density, s, spin kinetic energy density, T, and the 
full spin-current pseudotensor, J, as 

E = J(frH{p,T,j,s,T,3;r). (18) 

The time-odd terms (j, s, T) vanish for static calculations of even-even nuclei, while they are 
present for odd mass nuclei, in cranking calculations, as well as in TDHF. The spin-current 
pseudotensor, J, is time-even and does not vanish for static calculations of even-even nuclei. 
It has been shown [36l - l4T| that the presence of these extra terms are necessary for preserving 
the Galilean invariance and make an appreciable contribution to the dissipative properties of 
the collision. Our TDHF program includes all of the appropriate combinations of time-odd 
terms in the Skyrme interaction. In addition, commonly a pairing force is added to incor¬ 
porate pairing interactions for nuclei. The implementation of pairing for time-dependent 
collisions is currently an unresolved problem although small amplitude implementations ex¬ 
ist P^44| . However, for reactions with relatively high excitation this is not expected to be 
a problem. 


3 Numerical methods 

In this section we discuss the numerical details of performing TDHF calculations of nu¬ 
clear collisions as well as the density-constraint method, which is crucial for ion-ion potential 
calculations. Relatively recently it has become feasible, for the first time, to perform TDHF 
calculations on a 3D Cartesian grid without any symmetry restrictions and with much more 
accurate numerical methods |39l |45l |46]. At the same time the quality of effective interac¬ 
tions has also been substantially improved [47H50] . 


3.1 Discrete variation and lattice eqnations 

The lattice solution of differential equations on a discretized mesh of independent vari¬ 
ables may be viewed to proceed in two steps: (I) Obtain a discrete representation of the 
functions and operators on the lattice. (2) Solve the resulting lattice equations using itera¬ 
tive techniques. Step (I) is an interpolation problem for which we could take advantage of 
the techniques developed using the basis-spline functions [5T1|52]. The use of the basis-spline 
collocation method leads to a matrix-vector representation on the collocation lattice with a 
metric describing the transformation properties of the collocation lattice. 

In order to obtain a set of lattice equations which preserve the conservation laws associ¬ 
ated with the continuous equations it is essential to develop a modified variational approach. 
This goal is achieved by performing a variation to the discretized form of a conserved quan¬ 
tity, i.e. total energy. Consequently, the resulting equations will preserve all of the conserved 
quantities on the lattice. For the TDHF equations we consider a general discretized form of 
the action 


5 = 


dt 

\ 




dt 


(a/37) 



(19) 


where indices a, /3 ,7 denote the lattice points in three-dimensional space, and AVap^ is the 
corresponding infinitesimal volume element. Due to the presence of derivative operators in 
the Hamiltonian the explicit form of these expressions will depend non-locally on the lattice 
indices. The general variation, which preserves the properties of the continuous variation, 
is given by 


HljaPl) _ I 

SipKa'P'j') AI4/37 




( 20 ) 


Until recently, most HF and TDHF calculations have been performed using finite-difference 
lattice techniques. The details of the discrete variation for the finite-difference case are given 
in Refs. [351 [S3]. Below we outline a procedure for using the BSCM for the numerical solution 
of HF and TDHF equations. Further details of the BSCM is published elsewhere m- 






3.2 Basis-splines 

Given a set of points or knots denoted by the set {xi} a basis-spline (B-spline denoted by 
B^) function of order M is constructed from continuous piecewise polynomials of order M — 
1 [54] . B-splines have continuous derivatives up to (M — 2)”^* derivative and a discontinuous 
(M — 1)®* derivative. We only consider odd order splines or even order polynomials for 
reasons related to the choice of the collocation points. The B-spline is nonzero only 
in the interval (xi,Xi+M)- This property is commonly referred to as limited support. The 
knots are the points where polynomials that make up the B-spline join. In the interval 
containing the tail region B-splines fall off very rapidly to zero. An example of order M = 5 
splines extending over a physical region is illustrated in Figj^ We can also construct exact 
derivatives of B-splines provided the derivative order does not exceed M — 1. 



Figure 1: A region of space with physical boundaries located at knots xm and xm+n for M = 5 and 
N = 5. The B-spline which begins at the first knot xi has its tail in the physical region. The last 
B-spline which begins within the physical boundaries is It extends up to the last knot 3;jv+2A'/-i- 

A continuous function /(x), defined in the interval {xmin,Xmax), can be expanded in 
terms of B-spline functions as 


/(x)=^Bf(xK, 

i 


( 21 ) 


where quantities c* denote the expansion coefficients. We can solve for the expansion coef¬ 
ficients in terms of a given (or to be determined) set of function values evaluated at a set 
of data points more commonly known as collocation points. There are a number of ways 
to choose collocation points imiHj, however, for odd order B-splines a simple choice is to 
place one collocation point at the center of each knot interval within the physical boundaries 


Xa = 




a = 1,..., iV. 


( 22 ) 


Here, Xm = xn+m = x^ax, and N is the number of collocation points. We can now 

(23) 


write a linear system of equations by evaluating (21) at these collocation points 

fa='^ BaiC’- , 


where /„ = f{xa), and Bai = B^{xa). In order to solve for the expansion coefficients the 
matrix B needs to be inverted. However, as it stands matrix B is not a square matrix, since 
the total number of B-splines with a nonzero extension in the physical region is N + M — 

In order to perform the inversion we need to introduce additional linear equations which 
represent the boundary conditions imposed on f(x) at the two boundary points, xm and 
xm+n- The essence of the lattice method is to eliminate the expansion coefficients c® using 
this inverse matrix. The details of using the boundary conditions (or periodic boundaries) 










and inverting the resulting square matrix are discussed in Ref. m- Following the inversion 
the coefficients are given by 

c* = ^[B-Y“/a- (24) 

Ct 

One can trivially show that all local functions will have a local representation in the finite 
dimensional collocation space 

fix) U ■ (25) 

The collocation representation of the operators can be obtained by considering the action 
of an operator O onto a function /(a;) 

Of{x)=Y^[OB^{x)]c\ (26) 

i 

If we evaluate the above expression at the collocation points x^ we can write 

[OfU^Y.[OBUB . (27) 


Substituting from Eq. ( (24| for the coefficients B we obtain 

[o/]„ = 

ij3 

= Y.Oifp, ( 28 ) 

p 


where we have defined the collocation space matrix representation of the operator O by 


oi = Y.[obu [b- 


12/3 


(29) 


Notice that the construction of the collocation space operators can be performed once and 
for all at the beginning of a calculation, using only the given knot sequence and collocation 
points. Due to the presence of the inverse in Eq. (29) the matrix O is not sparse. In 


practice, operator O is chosen to be a differential operator such as d/dx or df /dx^. By a 
similar construction it is also possible to obtain the appropriate integration weights on the 
collocation lattice m- 


3.3 Discrete HF equations 


Since the detailed derivation of the BSCM representation of the TDHF equations involve 
many terms that are present in the energy functional, we will only show few terms. The 


three-dimensional expansion in terms of B-splines is a simple generalization of Eq. (21) 


'ipxix,y,z) = Biix)Bjiy)Bkiz) 


(30) 


ijk 


The knots and collocation points for each coordinate can be different. With the appropriate 
definition of boundary conditions all of the discretization techniques discussed in the previous 
section can be generalized to the three-dimensional space. The details of this procedure are 
given in [5T] . 

As an example for a local term let us consider a part of the to contribution to the total 
energy 


^(1 + y ) [ = ^(1 -p )^) ^ w°‘wf^w^[pial3-f)] 


(31) 


where on the right-hand side we have written the discretized form on a collocation lattice 
with collocation weights denoted by w. Here, a, j5 ,7 represent the collocation points in a;, j/. 




and z directions, respectively. In order to be able to perform the variation with respect to 
the single-particle states we rewrite equation explicitly 

a/37 


Using Eq.(20) in the variation of Eq.(32) we obtain (after replacing the primed indices with 
unprimed ones) the contribution 


to(l + y)p(a/37)V'A(a/37) , 


(33) 


where we have rewritten a summation as the total density. The same procedure can be 
carried out for the nonlocal terms in the energy density. A typical term is illustrated below 


(^V'a )a/37 


+ V'J (a/3'7) i 

a' /3' 




where the matrices D denote the first derivative matrices in x,y,z directions (they can be 
different although the notation does not make this obvious), calculated as described in the 
previous section. Finally, the HE equations can be written as matrix-vector equations on 
the collocation lattice 

hip^ —^ h • -0A • (34) 

The essence of this construction is that the terms in the single-particle Hamiltonian h are 
matrices in one coordinate and diagonal in others. Therefore, h need not be stored as a full 
matrix, which allows the handling of very large systems directly in memory. The details of 
this procedure are discussed below. 


3.4 Solution of the discrete equations 

In this subsection we will outline some of the numerical methods developed for the 
solution of the discretized HE and TDHF equations. The subsection is divided into two 
parts; the first part discusses the static iteration methods and the solution of the field 
equations and in particular the powerful damped relaxation method [55]. In addition we 
discuss the implementation of external constraints on the HE equations. The second part 
of the subsection introduces a number of time-evolution methods used in our calculations. 
Typical numerical accuracies are also discussed. 


3.4.1 Static solutions 


The solution of the HE equations (34) represent the problem of finding the few lowest 
eigenvalues of a very large Hamiltonian matrix. Furthermore, due to the fact that we are 
dealing with a self-consistent problem the matrix elements must be recalculated at every 
iteration. However, in practice the matrix elements need not be stored. Instead, one can 
make use of the inherent sparsity to dynamically construct the operation of the single-particle 
Hamiltonian onto a statevector. The basic operation is 


tp' = h ■ Ip 


(35) 


where the construction of the right hand side is done by explicitly programming the required 
linear combinations of the elements of ip to give xp'. In this approach the only storage 
requirements are for the statevectors and small matrices present in the Hamiltonian. 

The lattice equations are solved by using the damped relaxation method described in 
Refs. [291155] . A simple way for introducing the damped relaxation method is by pointing out 
its resemblance to the so-called imaginary time method. A more formal discussion is given 





in Ref. |55j where the generalization to the relativistic Dirac equation is also introduced. 
We start with the TDHF equations 


ih 


dt 


= ■ 


(36) 


In terms of the discretized time = nAt the solution at time n + 1 can be obtained from 
time n by (see below) 

^n+l ^ ^-^Athyh^n ^ ^ 37 ) 

where h" is the single-particle Hamiltonian at the iteration. The imaginary time-step 
method consists of the transformation At —?► —iAt 




(38) 


where xq = At/h, and we have taken out a trivial phase from ■0^. The expansion of the 
exponential to first order in xq yields the imaginary time iteration scheme 

Xl-^^=0[xl-xo{h--el)xl], (39) 

where O stands for Gram-Schmidt orthonormalization, which is necessary to ensure the 


orthonormality of the single-particle states at each iteration. In Eq. (39) the index n is 


no longer associated with time and it simply becomes an iteration counter. It is clear from 


Eq. (381 that the exponential acts as a filter in selecting the lowest eigenvalues of h and leads 


to the minimization of the HE energy. The generalization of the imaginary time method, 
where we introduce the damping matrix D, results in the damped relaxation method 


Xr = 0 [Xa - ^oD(Eo)(h" - eM ■ 


(40) 


The damping operator D is chosen to be 

d(e;o) = + ^ 

Eq 


r 

-1 


-1 


^o_ 


Eq 


O 


1 -1 


where T denotes the kinetic energy operator. Limits can be established for the ranges of the 
parameters xq and Eg |55j . but in practice fine-tuning is necessary for optimal performance. 
Two convergence criteria are used in practical calculations; one being the fractional change 
in the HE energy 

rnn+l rpn 

AE- = , (41) 


and other the fluctuations in energy 

77 = y< E2 > - < iJ >2 


(42) 


The fluctuations are a more stringent condition than the simple energy difference between 
two iterations. In practice, we have required rj to be less than 10“®. For this value of the 
energy fluctuation the fractional change in the HE energy is about 10“^^. 

The calculation of the HE Hamiltonian also requires the evaluation of the direct Coulomb 
contribution. However, since the calculation of the three-dimensional Coulomb integral is 
very costly, we solve instead the corresponding differential equation 

V^Ec(r) = -47re^pp(r) . (43) 

Details of solving the Poisson equation using the BSCM is given in Ref. [5T| . 














3.4.2 Constrained HF calculations 


It is sometimes desirable to solve the static HF equations away from the global minimum 
in energy. Such situations usually arise in the study of fission barriers and in the study of 
long-lived superdeformed states of nuclei that can be formed during low energy heavy-ion 
collisions. These methods have been instrumental for the understanding of the formation of 
nuclear molecules [53]. All of these cases require the existence of a stable minimum which 
does not coincide with the ground state configuration. The usual approach is to study the 
HF energy of a nuclear system by keeping certain macroscopic degrees of freedom at pre¬ 
specified values. This results in a multi-dimensional energy surface from which extremum 
values can be obtained. The reliability of these results depend strongly on the correct 
identification of the relevant macroscopic degrees of freedom. However, as we will see below 
a special constrained HF method, density constrained HF, has also been developed which 
allows the minimization of the energy along a TDHF trajectory. 

The goal is to devise an iteration scheme such that the expectation value of an arbitrary 
operator Q does not change from one static iteration to next 

E < >= E < XaIQIXa > • (44) 

A A 

Furthermore, we require this expectation value to be a fixed number Qq. A procedure can be 
developed by using Lagrange multipliers that are dynamically adjusted [29|. We start with 
the static HF iteration scheme modified by the addition of a constraint (we have omitted 
the damping matrix D in the equations below for simplicity) 

xl-^^ = 0[xl-xo{h- + XQ-el)xl]. (45) 

In Ref. [53] we give a set of exact equations which preserve the expectation of the constraining 
operator to order Xq. However, these equations involve the calculation of exchange terms 
and may become costly. Instead, one can develop a simpler iterative scheme as follows. 
Perform an intermediate step 

^n+l /2 ^ ^ ^ 

and calculate the difference 

5 Q-+ 1/2 = ^ < ^nH-l/ 2 |Q|^nH-l /2 > _ ^ > . ( 47 ) 

A A 

In analogy with the exact case the Lagrange parameter A is altered to reduce this difference 


\ n+l \ n I 

A ^ = A + Co - 




(48) 


2a;oEA < XaIQ^Ix” > 

where cq and do are empirical parameters replacing the exchange terms. In terms of these 
intermediate states the (n -I- l)st step is given by 


xT^ = o[x 


n-i- 1/2 


where 


<5A" = 


-xo(A”+'-A" + dA”)QxA 

Ea < xllQlxl > -Qo 


ti+ 1 / 2 i 


(49) 


(50) 


2a;oEA < XaIQ^Ix” > +do 

The extension of the method which allows the entire density to be constrained is straight¬ 
forward. In this case we would like to constrain a continuous density 


P"(r)=ElXA(r)|E 


(51) 


to be equal to / 0 o(i’ )• The constraining operator Q becomes the density operator p{r ) defined 
in the single-particle space 

<X^ld(r)|XA >= lx^(r)P 


( 52 ) 




and the product XQ is replaced by an integral 


XQ 


(frX{r)p{r) = A(r) . 


(53) 


The last equality is due to the fact that in coordinate space p(r) is a delta function. An 


iterative scheme for A"(r) is given by 


: n+l/2 

A"+i(r) = A"(r) + co„ ^ , , 

2xop^{r) + do 

(54) 

where 


dp"+i/2(r)=p"+i/2(r)-po(r) 

(55) 

is obtained from half-time iteration step 


^n+l/2 ^ ^ ^ ^ 

(56) 


Note that in 
iteration and 
written as 


these equations we require that the density remain equal to po (r) at every 
not just at the final step. Using these wavefunctions the full iteration can be 

- xo(A"+i(r) - A"(r) + 5A”(r))xr'/'] , (57) 


where 


5A"(r) = co 


p"(r) -po(r) 
2xoPo{y) + do 


(58) 


During the density constrained HF iterations the single-particle states readjust to mini¬ 
mize the energy while the initial density is kept fixed. In practical calculations the parameter 
xq has been replaced by the damping operator and the constants cq and do were chosen to 
be 1.9 and 5 x 10“®, respectively. 


3.4.3 Time evolution 

The formal solution of the TDHF equations ([^ is 

'^x{t) = Uit,to)'il^x{to) , (59) 


where we have omitted the spatial coordinates for simplicity and the time propagator U(t, to) 
is given by 


U{t,to) = T exp 



(60) 


The quantity T denotes time-ordering which is necessary in the general case. In practical 
calculations we discretize time as 


tn = nAt n = 0,1,2, ...N , (61) 

and express the evolution operator in successive infinitesimal pieces 

Uit,to) = U{t,tN-l)U{tN-l,tN-2)---U{ti,to) ■ (62) 

In this case the time-ordering operator can be ignored. For three-dimensional calculations 
the exponential operator is expanded as a Taylor series 

[/«„«, t„)» (^1+g ^ (63) 

The expansion of the operator requires repeated applications of h onto the wavefunctions. 
In practice, only 6 — 8 terms are needed for the conservation of the norm at 1 part in 10“^° 
level during the entire time-evolution. The expansion method is also attractive due to the 
fact that it only involves matrix vector operations which could be easily customized for 
vector or parallel computers. 







4 Numerical results 


4.1 Fusion Calculations for Neutron-Rich and Stable Nuclei 

In fusion experiments with neutron-rich radioactive ion beams, the dynamics of the 
neutron skin usually enhances the sub-barrier fusion cross section over that predicted by 
a simple static barrier penetration model, but in some cases suppression of fusion is also 
observed. Most recently, we have investigated sub-barrier fusion and pre-equilibrium giant 
resonance excitation between various tin + calcium isotopes [51 and calcium -|- calcium 
isotopes [20]. Finally, we have studied sub-barrier fusion reactions between both stable and 
neutron-rich isotopes of oxygen and carbon [58H60j that occur in the neutron star crust. In 
all cases, we have found good agreement between the measured fusion cross sections and the 
DC-TDHF results. This is rather remarkable given the fact that the only input in DC-TDHF 
is the Skyrme effective N-N interaction, and there are no adjustable parameters. 

In several experiments, large enhancements of sub-barrier fusion yields have been ob¬ 
served for systems with positive Q values for neutron transfer. Recently, a series of exper¬ 
iments has been carried out with radioactive ^^^Sn beams and with stable ^^^Sn beams on 
targets m- It turns out that the ^®Ca+Sn systems have many positive Q values for 
neutron-pickup while all the Q values for ^^CaH-Sn are negative. However, the data analysis 
reveals that the fusion enhancement is not proportional to the magnitudes of those Q values. 




Figure 2: Left: DC-TDHF calculation |56l I57| of heavy-ion interaction potentials for the neutron-rich 
system ^^^Sn-h‘^®Ca. The potentials are shown at five Ec.m. energies. Right: Total fusion cross section as a 
function of center-of-mass energy. The experimental data are taken from Ref. m- 

First, we consider the neutron-rich system ^^^Sn-t-'^^Ca. Results for the heavy-ion in¬ 
teraction potential V(R) with the DC-TDHF method [551 EZ] are shown in Fig. on the 
left side. In general, the ion-ion potentials are energy-dependent, and they are shown here 
at five Ec.m. energies. Our results demonstrate that in these heavy systems the potential 
barrier height increases dramatically with increasing energy, and the barrier peak moves 
inward towards smaller i?-values. On the right side of Fig. [^we show the calculated fusion 
cross section as a function of E^.m.- In this case, we have used the theoretical cross sec¬ 
tions obtained with the energy-dependent DC-TDHF potentials [55]. We can see that our 
theoretical cross sections agree remarkably well with the experimental data. 

Particularly puzzling is the experimental observation of a sub-barrier fusion enhance¬ 
ment in the system ^^^Sn-|-'*°Ca as compared to more neutron-rich system ^^^Sn-|-'*®Ca. 
This is difficult to understand because the 8 additional neutrons in "‘^Ca should increase 
the attractive strong nuclear interaction and thus lower the fusion barrier, resulting in an 
enhanced sub-barrier fusion cross section. But the opposite is found experimentally. A 
coupled channel analysis [61] of the fusion data with phenomenological heavy-ion potentials 
yields cross sections that are one order of magnitude too small at sub-barrier energies, de¬ 
spite the fact that these ion-ion potentials contain many adjustable parameters. Using the 
microscopic DC-TDHF approach we are able to explain the measured sub-barrier fusion en- 

















Figure 3: Left: Comparison of DC-TDHF heavy-ion interaction potentials, calculated at -Ec.m. = 120 
MeV, for the systems ^^^Sn+'^'^Ca and ^^^Sn-h^^Ca. Right: Total fusion cross sections for ^^^Sn-h^'^Ca as 
a function of center-of-mass energy. The cross sections have been calculated for several energy-dependent 
potentials. The experimental data are taken from Ref. m- 


hancement m in terms of the narrower width of the ion-ion potential for ^^^Sn-|-^*^Ca, while 
the barrier heights and positions are approximately the same in both systems, see Fig. 
The energy-dependence of the heavy-ion interaction potentials is crucial for understanding 
the low-energy fusion data. 

As another example, we discuss our theoretical analysis of recent -|-‘*°’"^®Ca fusion 

data [62] which extend older fusion data to low sub-barrier energies. Comparison of the 
sub-barrier cross sections with those calculated using standard coupled-channel calculations 
suggested a hindrance of the fusion cross-sections at deep sub-barrier energies. The DC- 
TDHF results for the ion-ion potentials and fusion cross section [20] are shown in Fig. 

In contrast to the heavy systems discussed above, we find that the ion-ion potentials for 



Figure 4: DC-TDHF ion-ion potentials (left) and fusion cross-sections (right) for various isotopes of the 
Ca-I-Ca system [20|. The experimental data are taken from Ref. m- 

lighter systems depend only weakly on the energy, and the potential barrier corresponding 
to the lowest collision energy gives the best fit to the sub-barrier cross-sections since this 
is the one that allows for more rearrangements to take place and grows the inner part of 
the barrier. Considering the fact that historically the low-energy sub-barrier cross-sections 
of the "^°Ca-|-'^®Ca system have been the ones not reproduced well by the standard models, 
the DC-TDHF results are quite satisfactory, indicating that the dynamical evolution of the 
nuclear density in TDHF gives a good overall description of the collision process. 

Figure shows the excitation energy, E*{R), for the three systems studies here. The 
excitation energy was calculated for the same value of e = Ec.m./h- = 2.75 MeV for all 
systems, which corresponds to collision energies of 55, 60, and 66 MeV, respectively. All 










curves initially behave in a similar manner, at large distances the excitation is zero, as 
the nuclei approach the barrier peak the excitations start and monotonically rise for larger 
overlaps. The interesting observation is that the excitations for the intermediary ^^Ca+'^^Ca 
system start at a slightly earlier time and rise above the other two systems. This may be 
largely due to the fact that an asymmetric system has some additional modes of excitation 
in comparison to the other two symmetric systems. 



Figure 5: Excitation energy E*{R), for three isotopic combinations of the Ca+Ca system 1201 . 


4.2 Hot and Cold Fusion to Produce Superheavy Elements 

In connection with superheavy element production, we have studied the hot fusion reac¬ 
tion ^®Ca-|-^^®U and the cold fusion reaction ^°Zn-|-^°®Pb (TH]. Considering hot fusion, "‘^Ca 
is a spherical nucleus whereas g large axial deformation. The deformation of 

strongly influences the interaction barrier for this system. This is shown in Fig. which 
shows the interaction barriers, V{R), calculated using the DC-TDHF method as a function 
of c.m. energy and for three different orientations of the nucleus. The alignment angle 
/3 is the angle between the symmetry axis of the nucleus and the collision axis. 



Figure 6: Left: Potential barriers, V{R), obtained from DC-TDHF calculations 1181 as a function of 
®c.m. energy and orientation angle f} of the ^SSpj nucleus. Also shown are the experimental c.m. energies. 
Right; Capture cross-sections as a function of Ec.m. energy (black circles). Also shown are the experimental 
cross-sections 16311641 (red squares). 

















The barriers for the polar orientation {/3 = 0°) of the nucleus are much lower 
and peak at larger ion-ion separation distance R. On the other hand, the barriers for the 
equatorial orientation (/3 = 90°) are higher and peak at smaller R values. We observe that at 
lower energies the polar orientation results in sticking of the two nuclei, while the equatorial 
orientation results in a deep-inelastic collision. We have also calculated the excitation energy 
E*{R) as a function of c.m. energy and orientation angle f3 of the nucleus. The system 
is excited much earlier during the collision process for the polar orientation and has a higher 
excitation than the corresponding collision for the equatorial orientation. 

To obtain the capture cross-section, we calculate potential barriers V{R,fi) for a set of 
initial orientations /3 of the nucleus. Then we determine partial cross sections cr(/3) and 
perform an angle-average, including the dynamic alignment arising from Coulomb excitation 
of the g.s. rotational band in the entrance channel. In Fig. we show our results for the 
capture cross-sections which are in remarkably good agreement with experimental data. 

One of the major questions that is asked by the experimental superheavy element com¬ 
munity is why a ^®Ca beam is so crucial in forming such systems and whether one could 
produce new super heavy nuclei using projectiles different than "^^Ca and actinide targets. 
Some possible projectiles include ®°Ti, ^'^Cr, ^®Fe and ®"^Ni. Several reactions with these 
projectiles were tried, no SHE event reported so far, only cross section limits. Typical re¬ 
actions studied already are 238pj (SHIP GSI), ^^®Cm-|-®^Cr (SHIP) and ^“^^Cm-l-^^Fe 

(Dubna), ^OTi-b^^SBk and soTi-b^^^Cf (TASCA GSI). During the ICFN5 meeting at Sanibel 
Island, Florida a particularly low cross section limit, of about 50 fb, was reported for the 
production of isotopes of new element 119 in the reaction between ®°Ti and ^^®Bk. The 
reaction ^®Ca -|-^‘^®Bk makes superheavy isotopes of element 117 with cross-sections of 2 — 3 
picobarns, by contrast the ®°Ti -|-^'*®Bk reaction yields a cross-section limit of 50 fb, about 
50 times lower. TDHF can be used to calculate capture cross-sections, excitation energies, 
and other dynamical properties for some of these systems in order to find the discerning 
physical property among them. 

4.3 Nuclear Astrophysics: Fusion of Oxygen and Carbon in the 
Neutron Star Crust 

Fusion of very neutron rich nuclei may be important to determine the composition and 
heating of the crust of accreting neutron stars. We have studied sub-barrier fusion reactions 
between both stable and neutron-rich isotopes of oxygen and carbon. In Fig. (left side) 
we show the DC-TDHF potential barriers for the C-bO system [5S]. The higher barrier 
corresponds to the ^^C-b system and has a peak energy of 7.77 MeV. The barrier for 
the ^^C-b system occurs at a slightly larger R value with a barrier peak of 6.64 MeV. 
The right side of Fig. shows the corresponding cross sections for the two reactions. Also 
shown are the experimental data from Ref. [SS]. The DC-TDHF potential reproduces the 
experimental cross-sections quite well for the ^^C-b system, and the cross section for the 
neutron rich ^^C-b^^O is predicted to be larger than that for ^^C -b^®0. Similar agreement 
with fusion data for the ^®0-b system and the predicted cross-sections for heavier oxygen 
isotopes were also calculated [55] . 


Summary 

In this manuscript we have outlined the microscopic study of fusion using the DC-TDHF 
approach. For a variety of heavy-ion reactions, we discuss the results of our numerical 
calculations of fusion / capture cross-sections at energies below and above the Coulomb 
barrier and compare these to experimental data if available. In this context we discuss a 
variety of stable and neutron-rich reaction partners, and we examine hot fusion reactions 
leading to the formation of superheavy elements. Fusion of very neutron rich light nuclei also 
plays a major role in nuclear astrophysics where it determines the composition and heating 
of the crust of accreting neutron stars. We have shown that microscopically obtained ion-ion 
potentials do give a reasonably good description of these fusion cross-sections. 




Figure 7: Left: DC-TDHF heavy-ion potentials for the systems Right: Total fusion cross 

sections versus center-of-mass energy for fusion of carbon with oxygen isotopes m- The experimental data 
are from Ref. m- 


The fully microscopic TDHF theory has shown itself to be rich in nuclear phenomena and 
continues to stimulate our understanding of nuclear dynamics. The time-dependent mean- 
field studies seem to show that the dynamic evolution builds up correlations that are not 
present in the static theory. While modern Skyrme forces provide a much better description 
of static nuclear properties in comparison to the earlier parametrizations there is a need 
to obtain even better parametrizations that incorporate deformation and reaction data into 
the ht process. 
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